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The effect of force field type on the predicted elastic properties of a polyimide is examined 
using a multiscale modeling technique. Molecular Dynamics simulations are used to predict 
the atomic structure and elastic properties of the polymer by subjecting a representative 
volume element of the material to bulk and shear finite deformations. The elastic properties 
of the polyimide are determined using three force fields: AMBER, OPLS-AA, and MM3. 
The predicted values of Young’s modulus and shear modulus of the polyimide are compared 
with experimental values. The results indicate that the mechanical properties of the 
polyimide predicted with the OPLS-AA force field most closely matched those from 
experiment. The results also indicate that while the complexity of the force field does not 
have a significant effect on the accuracy of predicted properties, small differences in the 
force constants and the functional form of individual terms in the force fields determine the 
accuracy of the force field in predicting the elastic properties of the polyimide. 


I. Introduction 

Nanostructured materials have the potential to provide significant increases in stiffness-to-weight and strength- 
to-weight ratios relative to current materials used for aerospace structural applications. To facilitate the 
development of these materials for this purpose, constitutive relationships must be developed that predict the bulk 
mechanical properties of the materials as a function of the molecular structure. 

One promising approach to accomplish this type of modeling over multiple length scales is to establish the use 
of molecular dynamics (MD) simulations to accurately predict performance at the atomistic level. MD simulations 
can be used to predict the equilibrated molecular structure of a material and the behavior of the molecular system 
when subjected to applied mechanical deformations. Many studies have focused on the modeling and simulation of 
polymers and polymer-based nanocomposites via MD techniques. 1 " 14 These studies have demonstrated that MD 
techniques can be effectively used to predict both structure and properties of polymer-based material systems. An 
important factor in the accurate prediction of properties of material systems via MD is the selection of the atomic 
potential. Several simplified atomic potentials, or force fields, for organic-based systems have been developed in 
recent years that describe the interactions between bonded and non-bonded atoms. 15 ' 36 Each of these force fields has 
been characterized via experimental techniques and quantum computations and has their own set of unique 
parameters and functional forms. Even though it is expected that these different parameters and forms will affect the 
relationship between force field type and predicted mechanical properties, little is known about this cause-and-effect 
relationship. 

The objective of the present paper is to use MD simulations with different force fields to predict elastic 
properties of a polymer system. The predicted values of elastic properties from each force field will be compared to 
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experimentally measured values. The polymer system is a polyimide from 3,3’,4,4’-biphenyltetracarboxylic 
dianhydride (BPDA) and l,3-bis(4-aminophenoxy)benzene (APB) monomers (Figure 1). 37,38 The three force fields 
used are the AMBER, OPLS-AA, and MM3 force fields. From these comparisons, the most appropriate force field 
for the prediction of mechanical properties of polymer-based nanocomposite systems is determined. 



Figure 1. Schematic illustration of BPDA (1,3,4) APB polyimide monomer unit 


II. Force Fields 


Three force fields were used in this study to simulate the polymer molecule deformation and compute the 
corresponding mechanical properties; the AMBER 18,24 ' 27 (without electrostatic interactions), OPLS-AA, 28,29 and 
MM3 15,30 " 36 force fields. Each of the force fields has a unique functional form and set of force constants. The total 
potential energy of a simulated molecular system computed with the AMBER force field is 
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where the summations are taken over all of the corresponding interactions in the molecular model; K r and Kg are the 
bond-stretching and bond-angle bending force constants, respectively; r and r eq are the bond length and equilibrium 
bond length, respectively; 9 and 9 eq are the bond angle and equilibrium bond angle, respectively; V n /2 , C n, and (j) 
are the torsion magnitude, phase offset, periodicity of the torsion, and the torsion angle, respectively; and £ u , r u , and 
a jj are van der Waals well depth, non-bonded distance between atoms / and J, and the equilibrium distance between 
atoms / and J , respectively. 

The total potential energy of the molecular model computed with the OPLS-AA force field is 
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where qi is the partial charge of atom /, e is the elementary charge, and the remaining quantities are those already 
defined for the AMBER force field. 

For the MM3 force field, the total potential energy is 
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where K r0 , K i/:r , and K m , are force constants; r' and r' are the bond length and equilibrium bond length, 
respectively, of the adjacent covalent bond; and O' and O' are the bond angle and equilibrium bond angle, 

respectively, of the adjacent bond angle. The energy contribution from electrostatic forces, A e i e ctrostatic, is determined 
by either partial charges or dipole moments. The energies associated with all remaining non-bonded interactions, 
such as hydrogen bonding, is incorporated in A nb . Each of the force constants for these three force fields is unique 
for each force field and interacting atom types. For all three force fields, it was assumed that torsional force 
constants that were not defined in the respective literature references or by the simulation software 39 were zero- 
valued. 

TIT. Equivalent-Continuum Modeling 

The nonlinear-elastic (hyperelastic) properties of the four material systems were determined using the 
Equivalent-Continuum Modeling method. 8 ' 9 ' 40 41 This approach consisted of three steps. First, a representative 
volume elements (RVEs) of the molecular structure for each force field were chosen that accurately described the 
bulk structure of the material. Next, a constitutive law that described the behavior of the equivalent-continuum 
model was established. Finally, the energies of deformation of the two models were equated under identical sets of 
boundary conditions to determine each of the material parameters in the constitutive equation. Each of these steps is 
described in detail below. 

A. Representative Volume Element 

The RVEs of the molecular models were the equilibrium molecular structures of the polyimide for each force 
field determined using MD simulations. The RVE geometry of the molecular model selected was a cubic box. For 
this study, the molecular model of the polyimide was prepared with the aid of a reverse-mapping procedure from a 
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coarse-grained model. 42 This method is similar to those previously employed for multiscale polymer modeling. 43 ' 44 
A brief summary of the method is given here with additional details presented elsewhere. 46 

For each polyimide molecule, a linked vector model was used to represent the rigid rings that comprise the 
polyimide backbone (Figure 2). The linked vectors followed the contour of the molecule. The parameters used for 
this model consisted of angular distributions between consecutive vectors and long-range forces between beads 
placed along the midpoint of each vector. These parameters were estimated from molecular dynamics (MD) 
simulation of the polyimide monomers with the CVFF force field. 45 The centroids of the beads placed at the 
midpoint of each vector were the centers for interaction forces between non-adjacent beads along the chain of the 
polymer and between beads on different chains. 



Figure 2. Depiction of the mapping of the atomistic polymer model to the coarse-grained linked vector model 

Once the coarse-grained model for the single polyimide molecule was established, it was subsequently used to 
assemble the coarse-grained bulk model with multiple polymer molecules and the nanoparticle. The coarse-grained 
polymers were initially placed as random walk chains inside a simulation box (with periodic boundary conditions) 
close to their bulk density. In this initial placement, only the angular distributions between adjacent vectors along 
the chain were considered. Monte Carlo simulation was used to equilibrate the chains from their initial starting 
configuration. The bulk polymer model consisted of seven chains of polymers each composed of ten of the repeat 
units shown in Figure 1. This chain length is typical for MD simulations of polymers. The simulation box was a 
cube of side length 42 A. The periodic box dimensions were chosen to allow the polymer to be close to the 
equilibrium bulk density. The simulation ran at 650 K until relaxation of the autocorrelation function of the end 
vectors was achieved and the average centers of mass were displaced a distance greater than the average radii of 
gyration squared. After sufficient equilibration with the coarse-grained Monte Carlo model, the chains were 
reverse-mapped by replacing the deleted atoms back into position along the vectors of the coarse-grained model. 



Figure 3. RVE of the polyimide material 
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The resulting atomistic structures were subsequently minimized by the following procedure. A short energy 
minimization was applied to the structures. This was followed by an nPT (constant number of atoms, pressure, and 
temperature) MD simulation for 200 ps at 300 K and 1 atm using the AMBER, OPLS-AA, and MM3 force fields. 
This constant-pressure MD simulation allowed the atomistic structures to relax to the equilibrium density. The 
employed algorithm preserved the cubic structure of the simulation box while allowing the simulation box size to 
change. The final periodic boundary box sizes varied from 37.7 A to 47.2 A on a side depending on the force field 
used. After the nPTMD simulations, the densities of the pure polyimide matrix was about 1.3 g/cm 3 . This density 
is a reasonable value for polyimides. 47 An example of a RYE of the polyimide is shown in Figure 3. 


B. Constitutive Equation 

For the computational simulation of polymer materials subjected to finite deformations, it is assumed that the 
strain energy function is associated with stress and strain tensors that are thermodynamic work conjugates in the 
balance of mechanical energy. In particular, the second Piola-Kirchhoff stress tensor is 


S = 


(E) 

3E 


(15) 


where E is the Lagrangian strain tensor and V P C is the scalar strain-energy function of the equivalent continuum 
(defined per unit volume of the RVE in the reference configuration). The second Piola-Kirchhoff stress and 
Lagrangian strain tensors are henceforth referred to as the stress and strain tensors, respectively. Perhaps the 
simplest form of the strain energy function for isotropic hyperelastic materials is given by the classical Saint-Venant 
Kirchhoff model 


'P c (E) = |(trE) 2 + / /tr(E 2 ) (16) 

where A and // are the usual Lame constants, trE is the trace of the tensor E, and tr(E 2 ) is the trace of tensor E 2 . It 
can be shown 48 
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where I is the second-order identity tensor. Substitution of Equations (16) and (17) into Equation (15) yields 

S = Ttr(E)I + 2//E (18) 

Equation (18) is the stress-strain relationship of the isotropic polymer material. This relationship describes the 
mechanical behavior of the equivalent-continuum model. At this point, however, the materials constants A and ft are 
unknown, and must be determined using on the molecular structure of the polymer. This is accomplished in the 
proceeding modeling step. 


C. Energy Equivalence 

The energies of deformation of the equivalent-continuum, 'Pc, and molecular models, V P,„, were equated for 
identical sets of boundary conditions to determine the mechanical properties of the polyimide for each of the force 
fields. The molecular energy of deformation (potential energy) is 

, P„4(A'--A-) (19) 
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where A° lolal and A lotal are the potential energies of the molecular model from Equations (1), (2), and (7) before and 
after deformation, respectively, and V 0 is the initial volume of the RVE. For finite deformations, the kinematic 
equations of motion that describe the deformation of the boundary of the RVE are generalized by 

x = x(X,t) (20) 

where X are the material coordinates (before deformation, defined in Figure 3), x are the spatial coordinates (after 
deformation), and t is the time. The resulting strain that results from the deformation in Equation (20) is 

E = i(F r F-l) (21) 

where the superscript T denotes a tensor transpose and F is the deformation gradient tensor whose components are 
given by 



( 22 ) 


Because the polymer materials are assumed to be isotropic, only two independent material parameters are required 
to completely describe the three-dimensional linear-elastic constitutive behavior described by Equation (18). A 
single set of deformations is required to determine each of the two independent parameters. The two independent 
parameters that were determined for these polymers where the bulk modulus, /c and shear modulus, //. States of 
pure dilatation and pure shear were applied to the simulation box to determine k and //, respectively. 

For the bulk deformations, the kinematic equations of motion (Equation (20)) become 
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where cq, a 2 , cq, and a 4 are scalar constants and the superscripts indicate the deformation step. The coordinates x (1) , 
x <2) , x (3) , and x <4) correspond to bulk strains of 0.25%, 0.50%, 0.75%, and 1.0%, respectively. The relative 
deformation gradients (indicated by the apostrophe), which relate the deformation at a given strain level to those of 
the previous strain level are 
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Therefore, the deformation gradients that relate the coordinate for each strain level to those of the material 
coordinate system are 
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Using Equations (22) - (25), the constants oci, a 2 , a 3 , and a 4 where adjusted to achieve the exact desired strain levels 
of 0.25%, 0.5%, 0.75%, and 1.0% in Equation (21). These values are listed in Table I. 


Table I. Values of deformation parameters 


Deformation 

parameter 

Value (unitless) 

Finite-valued strain 
components 

a\ 

1.002497 

E u = E 22 = E 33 = 0.25% 

a 2 

1.002485 

E\\ = E 2 2 = £"33 = 0.50% 

a 2 

1.002472 

E\\ = E 22 = £"33 = 0.75% 

a 4 

1.002460 

E n = E 2 2 = E 33 =\.00% 

A 

0.001249 

Yu = Yn = Y \2 = 0.25% 

A 

0.001246 

Yu = Yn = Y\2 = 0.50% 

A 

0.001243 

Y 22 = Yn = Y \2 = 0.75% 

A 

0.001240 

Y 22 = 713 = Yn = 1-00% 


For the pure shear simulations, the deformations of the four strain levels are 
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where /3\, /T, /T, and fk are scalar constants and the superscripts 1, 2, 3, and 4 (k = 2, 3, and 4) correspond to pure 
shear strain levels of y 2 i = Yn = Yu = 0.25%, 0.50%, 0.75%, and 1.0%, respectively (/y = 2 Ey when i ± ./)• Similar to 
the case of the pure bulk deformation, the constants ft i, /%, Pi, and /? 4 were adjusted such that these shear strains 
were achieved by using Equations (21) and (24) - (26). The values of the /? constants are listed in Table 1. 

The energies of deformation of the molecular models were determined using MD simulations. Both bulk and 
shear deformations were applied to the equilibrium molecular structures for each force field by deforming the MD 
simulation boxes and all of the atoms in the models according to the applied deformation field. An nVT (constant 
number of atoms, volume, and temperature) simulation with periodic boundary conditions was subsequently used 


7 

American Institute of Aeronautics and Astronautics 




for each deformation to allow the box dimensions to remain fixed while the atoms were allowed to move into new 
equilibrium positions. The simulations were run up to 40 ps at 298 K, and were performed using the TINKER 
modeling package. 39 The potential energies of deformation of the molecular models were averaged over the final 10 
ps of each simulation, and the standard error of the energy fluctuations was recorded. A linear least-squares 
regression was performed for the potential energy of deformation over all four applied strain levels for each loading 
condition and force field. The elastic constants were subsequently determined using this data and Equation (16). 

For the MD simulations, periodic boundary conditions were applied such that atoms were free to cross the 
boundary of the deformed and undeformed simulation cells. Atoms that crossed the boundary entered the simulation 
cells on the opposite side, as described in detail elsewhere. 49 Therefore, none of the atoms in the MD simulations 
were kinematically over-constrained, as can occur in simulations of RVEs of heterogeneous material systems. 50 

IV. Results and Discussion 

Simulated values of the Young’s and shear moduli of the polyimide are shown in Tables 2 and 3, respectively, 
for the three force fields. As a reference, an experimentally-obtained value of Young’s modulus is shown in Table 
2. 37 The experimental shear modulus shown in Table 3 was obtained from the experimental Young’s modulus and 
an assumed Poisson’s ratio of 0.4. For each simulated value of Young’s and shear modulus, the percentage 
difference with respect the experimental value is also shown in Tables 2 and 3, respectively. 


Table II. Predicted and Experimental Young’s Modulus of Polyimide 


Method 

Young’s Modulus (GPa) 

% difference from 
experiment 

Experiment 

3.6* 

- 

Simulation 

12.7 

253% 

(AMBER force field) 

Simulation 

(OPLS-AA force field) 

2.7 

25% 

Simulation 
(MM3 force field) 

5.9 

64% 


* From reference [37] 


From Table 2, it is clear that the force field which yielded a simulated Young’s modulus that was closest to the 
experimental value was the OPLS-AA force field, with a value that was less than the experimental value. The 
simulations with the MM3 force field predicted a Young’s modulus that was larger than the experimental value. 
The AMBER force field yielded a Young’s modulus that significantly overestimated the experimental Young’s 
modulus. From Table 3, the same trends exist for the simulated shear moduli. Again, The OPLS-AA force field 
predicted a shear modulus below the experimental value, whereas the shear modulus associated with the MM3 force 
field is larger than the experimental value. The shear modulus from the simulations with the AMBER force field 
significantly overestimated the experimental value. 

Another important point can be drawn from the data in Tables 2 and 3. The functional forms of the AMBER and 
OPLS-AA force fields are very similar, as is evident from Eqns. (l)-(6). In addition, the force constants in these two 
force fields are also very similar. The functional forms of the AMBER and OPLS-AA are significantly more simple 
than that of the MM3 force field, Eqns (7)-( 14). The MM3 force field contains many nonlinear and cross-interaction 
terms. Because of the similarity in the results of the OPLS-AA and MM3 force field simulations, and the significant 
overestimations of the AMBER simulations, it appears that small differences in functional form of the force field or 
values of the force constants have a significant effect on simulated elastic properties of this polyimide. The 
complexity of the MM3 force field doesn’t appear to add significant accuracy for the specific simulations in this 
study. 
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Table III. Predicted and Experimental Shear Modulus of Polyimide 


Method 

Shear Modulus (GPa) 

% difference from 
experiment 

Experiment 

1.3* 

- 

Simulation 

6.0 

362% 

(AMBER force field) 

Simulation 

(OPLS-AA force field) 

0.9 

31% 

Simulation 
(MM3 force field) 

2.1 

62% 


* From reference [37] assuming v= 0.4 


The results of the AMBER and MM3 simulations in this study are similar to those of other studies in which the 
predicted values of mechanical properties are generally larger than the experimentally-determined properties. 13 ' 14 
One explanation of this trend is that RVEs in molecular dynamics simulations represent a nearly perfect molecular 
structure. In the actual experimental test specimens, the material contains low volume fractions of air pockets, 
inclusions, and unreacted monomer. Therefore, from this perspective, it is expected that simulated mechanical 
properties should be larger than the experimental properties. Even though the OPLS-AA force field predicted 
mechanical properties that were closest to experimentally-determined properties, it may underestimate the 
mechanical properties of the specific RVE of the polyimide material that was modeled. 

V, Summary 

In this study, the effect of force field type on the predicted elastic properties of a polyimide has been examined. 
An equivalent-continuum modeling approach was used to predict the elastic properties based only on the atomic 
structure of the polyimide. Molecular Dynamics simulations were used to predict the atomic structure of the 
polymer and the corresponding potential energy changes when subjected to applied finite deformations. The energy 
changes were equated with the strain energy of the equivalent-continuum model, which was determined using a 
thermodynamic potential function and a finite deformation framework. As a result, the elastic properties of the 
polyimide were determined using three force fields: AMBER, OPLS-AA, and MM3. The predicted values of 
Young’s modulus and shear modulus were compared with experimental values for the same polyimide. 

The results of this study indicate that the mechanical properties of the polyimide predicted with the OPLS-AA 
force field most closely matched those from experiment. While the predicted properties using the MM3 force field 
were also reasonably close to experimental properties, those of the AMBER force field significantly overestimated 
the mechanical properties of the polymer. These results also indicate that the complexity of the force field does not 
have a significant effect on the accuracy of predicted properties. Apparently, small differences in the force constants 
and the functional form of individual terms in the force fields determine the accuracy of the force field in predicting 
the elastic properties of the polyimide. The results of the simulations with the MM3 and AMBER force fields 
exhibit similar trends those of other studies. For these force fields, the mechanical properties determined through 
the modeling overestimate those from experiment. This is most likely due to the imperfections that are present in 
bulk polymer specimens which are tested experimentally, which cannot be easily modeled on the molecular scale. 
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